24. 3D Solid Mechanics

The equations of elasticity describe the deformation of solids due to applied forces.

from netgen.occ import *
from ngsolve.webgui import Draw
import ngsolve

box = Box((0,0,0), (3,0.6,1))
box.faces.name="outer"
cyl = sum( [Cylinder((0.5+i,0,0.5), Y, 0.25,0.8) for i in range(3)] )
cyl.faces.name="cyl"
geo = box-cyl

Draw(geo);

find edges between box and cylinder, and build chamfers (requires OCC 7.4 or newer):

cylboxedges = geo.faces["outer"].edges * geo.faces["cyl"].edges
cylboxedges.name = "cylbox"
geo = geo.MakeChamfer(cylboxedges, 0.03)

name faces for boundary conditions:

geo.faces.Min(X).name = "fix"
geo.faces.Max(X).name = "force"

Draw(geo);

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.1)).Curve(3)
Draw (mesh);

24.1. Linear elasticity

Displacement: \(u : \Omega \rightarrow R^3\)

Linear strain: $\( \varepsilon(u) := \tfrac{1}{2} ( \nabla u + (\nabla u)^T ) \)$

Stress by Hooke’s law: $\( \sigma = 2 \mu \varepsilon + \lambda \operatorname{tr} \varepsilon I \)$

Equilibrium of forces: $\( \operatorname{div} \sigma = f \)$

Displacement boundary conditions: $\( u = u_D \qquad \text{on} \, \Gamma_D \)$

Traction boundary conditions: $\( \sigma n = g \qquad \text{on} \, \Gamma_N \)$

24.2. Variational formulation:

Find: \(u \in H^1(\Omega)^3\) such that \(u = u_D\) on \(\Gamma_D\) $\( \int_\Omega \sigma(\varepsilon(u)) : \varepsilon(v) \, dx = \int_\Omega f v dx + \int_{\Gamma_N} g v ds \)\( holds for all \)v = 0\( on \)\Gamma_D$.

E, nu = 210, 0.2
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))

def Stress(strain):
    return 2*mu*strain + lam*Trace(strain)*Id(3)    

fes = VectorH1(mesh, order=3, dirichlet="fix")
u,v = fes.TnT()
gfu = GridFunction(fes)

with TaskManager():
    a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)
    pre = Preconditioner(a, "bddc")
    a.Assemble()

force = CF( (1e-3,0,0) )
f = LinearForm(force*v*ds("force")).Assemble()

from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates='\r', tol=1e-8)
gfu.vec.data = inv * f.vec

CG iteration 1, residual = 0.00018071481331955388     
CG iteration 2, residual = 7.786988796915474e-05     
CG iteration 3, residual = 8.498290164433323e-05     
CG iteration 4, residual = 6.815633040170596e-05     

CG iteration 5, residual = 6.11775382692514e-05     
CG iteration 6, residual = 4.649770504420967e-05     
CG iteration 7, residual = 3.310779730402245e-05     
CG iteration 8, residual = 2.7344520923959155e-05     

CG iteration 9, residual = 2.0649562682632237e-05     
CG iteration 10, residual = 1.6898694209669794e-05     
CG iteration 11, residual = 1.167617783759093e-05     
CG iteration 12, residual = 8.717921240061112e-06     

CG iteration 13, residual = 6.578387898741831e-06     
CG iteration 14, residual = 5.153749466780123e-06     
CG iteration 15, residual = 3.6889549037401625e-06     
CG iteration 16, residual = 3.0209233346364507e-06     

CG iteration 17, residual = 2.235745334989896e-06     
CG iteration 18, residual = 1.6042952744390234e-06     
CG iteration 19, residual = 1.1956102419098305e-06     
CG iteration 20, residual = 8.396604190595384e-07     

CG iteration 21, residual = 6.626002655213012e-07     
CG iteration 22, residual = 4.7361418022819315e-07     
CG iteration 23, residual = 3.6535329314711184e-07     
CG iteration 24, residual = 2.5398536815595273e-07     

CG iteration 25, residual = 2.0630254971094887e-07     
CG iteration 26, residual = 1.3934313082420504e-07     
CG iteration 27, residual = 1.11786406124239e-07     
CG iteration 28, residual = 7.469243140049882e-08     

CG iteration 29, residual = 5.5192257477400284e-08     
CG iteration 30, residual = 4.08068595089225e-08     
CG iteration 31, residual = 2.6842187937719415e-08     
CG iteration 32, residual = 2.2490719061315368e-08     

CG iteration 33, residual = 1.5197915668884992e-08     
CG iteration 34, residual = 1.0348799492401669e-08     
CG iteration 35, residual = 8.437759116603183e-09     
CG iteration 36, residual = 5.295435578134984e-09     

CG iteration 37, residual = 4.010668347750459e-09     
CG iteration 38, residual = 2.507359927456423e-09     
CG iteration 39, residual = 1.7886079437646787e-09     
CG iteration 40, residual = 1.4024020379671437e-09     

CG iteration 41, residual = 1.2104626173537914e-09     
CG iteration 42, residual = 9.819100290241733e-10     
CG iteration 43, residual = 5.789090852144978e-10     
CG iteration 44, residual = 4.0602644392311763e-10     

CG iteration 45, residual = 2.9437821092438174e-10     
CG iteration 46, residual = 2.4455861313904376e-10     
CG iteration 47, residual = 2.1605231661949959e-10     
CG iteration 48, residual = 1.3665589382993458e-10     

CG iteration 49, residual = 9.018142537816051e-11     
CG iteration 50, residual = 6.74360792970081e-11     
CG iteration 51, residual = 4.407374095862673e-11     
CG iteration 52, residual = 3.011196286661374e-11     

CG iteration 53, residual = 2.2418518465376623e-11     
CG iteration 54, residual = 2.1126709183253153e-11     
CG iteration 55, residual = 1.3136508225615745e-11     
CG iteration 56, residual = 9.137685996018394e-12     

CG iteration 57, residual = 6.600588526911087e-12     
CG iteration 58, residual = 4.418761686359022e-12     
CG iteration 59, residual = 2.9127641240192632e-12     
CG iteration 60, residual = 2.0525902224622983e-12     

CG iteration 61, residual = 1.3699881584818951e-12     
CG converged in 61 iterations to residual 1.3699881584818951e-12

with TaskManager():
    fesstress = MatrixValued(H1(mesh,order=3), symmetric=True)
    gfstress = GridFunction(fesstress)
    gfstress.Interpolate (Stress(Sym(Grad(gfu))))

Draw (5e4*gfu, mesh);

Draw (Norm(gfstress), mesh, deformation=1e4*gfu, draw_vol=False, order=3);